Immune population analysis

General Setup

Setup chunk

Setup reticulate

knitr::opts_chunk$set(fig.width = 8)
knitr::opts_knit$set(root.dir = normalizePath(".."))
knitr::opts_knit$get("root.dir")

[1] "/nas/groups/treutlein/USERS/tomasgomes/projects/liver_regen"

Load libraries

library(reticulate)
knitr::knit_engines$set(python = reticulate::eng_python)
py_available(initialize = FALSE)

[1] TRUE

use_python(Sys.which("python"))
py_config()

python:         /home/tpires/bin/miniconda3/bin/python
libpython:      /home/tpires/bin/miniconda3/lib/libpython3.8.so
pythonhome:     /home/tpires/bin/miniconda3:/home/tpires/bin/miniconda3
version:        3.8.3 (default, May 19 2020, 18:47:26)  [GCC 7.3.0]
numpy:          /home/tpires/bin/miniconda3/lib/python3.8/site-packages/numpy
numpy_version:  1.18.5

NOTE: Python version was forced by RETICULATE_PYTHON

Load data (from all cells)

library(Seurat)
library(ggplot2)
library(destiny)
library(plyr)
library(dplyr)
library(RColorBrewer)

Immune cell subsetting

Subset immune cells

allcells_css = readRDS(file = "data/processed/allcells_css.RDS")

Lymphoid cell analysis

Subset Lymphoid cells

DimPlot(all_imm_cells, reduction = "umap", group.by = "Condition")
DimPlot(all_imm_cells, reduction = "umap", group.by = "allcells_clusters")

DimPlot(all_imm_cells, reduction = "umap", group.by = "Donor")

DimPlot(all_imm_cells, reduction = "umap", group.by = "Phase")

FeaturePlot(all_imm_cells, reduction = "umap", features = c("MKI67", "ALB", "S100A8", "COLEC11",
                                                            "CLEC1B", "CLEC14A", "PTPRC", "EPCAM",
                                                            "CD14", "TRAC", "CD3E", "CD8A"))

Cluster and get markers

all_l_cells = RunUMAP(all_l_cells, dims = 1:25, verbose = F)
DimPlot(all_l_cells, reduction = "umap", group.by = "Condition")
DimPlot(all_l_cells, reduction = "umap", group.by = "allcells_clusters")

DimPlot(all_l_cells, reduction = "umap", group.by = "Donor")

DimPlot(all_l_cells, reduction = "umap", group.by = "Phase")

FeaturePlot(all_l_cells, reduction = "umap", features = c("MKI67", "ALB", "S100A8", "COLEC11",
                                                            "NKG7", "TRGC1", "PTPRC", "EPCAM",
                                                            "CD14", "TRAC", "CD3E", "CD8A"))

Add annotations

DimPlot(all_l_cells, reduction = "umap", group.by = "pca25_res.0.4", label = T)
DimPlot(all_l_cells, reduction = "umap", group.by = "Donor", label = F)

DimPlot(all_l_cells, reduction = "umap", group.by = "Condition", label = F)

Subset T/NK cells

new_l_labs = c("0" = "ab-T cells 1",
               "1" = "NK/gd-T cells",
               "2" = "ab-T cells 2",
               "3" = "Infiltrating NK cells", # PTGDS,CX3CR1 (infilt)
               "4" = "IgA+ Plasma cells",
               "5" = "B cells",
               "6" = "IgG+ Plasma cells",
               "7" = "Dividing NK cells",
               "8" = "ab-T cells (stress)")

all_l_cells$lymphoid_annot = new_l_labs[as.character(all_l_cells$pca25_res.0.4)]

Cluster and get markers

immune_pops = c("ab-T cells 1", "ab-T cells 2", "NK/gd-T cells", "Infiltrating NK cells")
all_t_cells = all_l_cells[,all_l_cells@meta.data$lymphoid_annot %in% immune_pops]
all_t_cells = suppressWarnings(SCTransform(all_t_cells, do.correct.umi = T, verbose = F, 
                                             vars.to.regress=c("unique_name","nCount_RNA"),
                                             variable.features.rv.th = 1, seed.use = 1,
                                             return.only.var.genes = F, 
                                             variable.features.n = NULL))

all_t_cells = RunPCA(all_t_cells, verbose = F)
all_t_cells = RunUMAP(all_t_cells, dims = 1:25, verbose = F)
DimPlot(all_t_cells, reduction = "umap", group.by = "Condition")

DimPlot(all_t_cells, reduction = "umap", group.by = "allcells_clusters")

DimPlot(all_t_cells, reduction = "umap", group.by = "lymphoid_annot")

DimPlot(all_t_cells, reduction = "umap", group.by = "Donor")

DimPlot(all_t_cells, reduction = "umap", group.by = "Phase")

FeaturePlot(all_t_cells, reduction = "umap", features = c("MKI67", "ALB", "S100A8", "COLEC11",
                                                            "NKG7", "TRGC1", "PTPRC", "EPCAM",
                                                            "CD14", "TRAC", "CD3E", "CD8A"))

# https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5007630/
new_t_labs = c("0" = "NK cells 1",
               "1" = "TRM cells", # may have CD4 and CD8; ITGA1
               "2" = "MAIT cells 2", # CXCR6, CCR6, CCR5, some RORC (not DE with the ILC3), some ZBTB16
               "3" = "NK cells 2",
               "4" = "NK cells 3",
               "5" = "CD8 ab-T cells 2",
               "6" = "MAIT cells 1", # hard to be sure, but has many hallmarks, even higher CD8A
               "7" = "Infiltrating NK cells",
               "8" = "Naive CD4+ T cells", 
               "9" = "gd-T cells", # CD3+ vs cl7, TRDC/TRG
               "10" = "CD8 ab-T cells 3",
               "11" = "CD8 ab-T cells 1",
               "12" = "Treg",
               "13" = "ILC3") # KIT, AHR, RORC, no CD3/CD8/CD4

all_t_cells$t_annot = new_t_labs[as.character(all_t_cells$pca25_res.0.9)]

Myeloid cell analysis

Subset Myeloid cells

# https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5007630/
new_t_labs = c("0" = "NK cells 1",
               "1" = "TRM cells", # may have CD4 and CD8; ITGA1
               "2" = "MAIT cells 2", # CXCR6, CCR6, CCR5, some RORC (not DE with the ILC3), some ZBTB16
               "3" = "NK cells 2",
               "4" = "NK cells 3",
               "5" = "CD8 ab-T cells 2",
               "6" = "MAIT cells 1", # hard to be sure, but has many hallmarks, even higher CD8A
               "7" = "Infiltrating NK cells",
               "8" = "Naive CD4+ T cells", 
               "9" = "gd-T cells", # CD3+ vs cl7, TRDC/TRG
               "10" = "CD8 ab-T cells 3",
               "11" = "CD8 ab-T cells 1",
               "12" = "Treg",
               "13" = "ILC3") # KIT, AHR, RORC, no CD3/CD8/CD4

Clustering and markers

immune_pops = c("cDCs 2", "Macrophages", "Kupffer cells", "cDCs 1", "pDCs")
all_m_cells = allcells_css[,allcells_css@meta.data$allcells_clusters %in% immune_pops]
all_m_cells = suppressWarnings(SCTransform(all_m_cells, do.correct.umi = T, verbose = F, 
                                             vars.to.regress=c("unique_name","nCount_RNA"),
                                             variable.features.rv.th = 1, seed.use = 1,
                                             return.only.var.genes = F, 
                                             variable.features.n = NULL))

all_m_cells = RunPCA(all_m_cells, verbose = F)
all_m_cells = RunUMAP(all_m_cells, dims = 1:25, verbose = F)
DimPlot(all_m_cells, reduction = "umap", group.by = "Condition")

DimPlot(all_m_cells, reduction = "umap", group.by = "allcells_clusters")

DimPlot(all_m_cells, reduction = "umap", group.by = "Donor")

DimPlot(all_m_cells, reduction = "umap", group.by = "Phase")

FeaturePlot(all_m_cells, reduction = "umap", features = c("MKI67", "ALB", "S100A8", "COLEC11",
                                                            "NKG7", "TRGC1", "PTPRC", "EPCAM",
                                                            "CD14", "TRAC", "CD3E", "CD8A"))

Add annotations

mrks_q = SoupX::quickMarkers(all_m_cells@assays$SCT@counts,
                             all_m_cells@active.ident, N = 10)
View(mrks_q[mrks_q$qval<=0.05,])

new_m_labs = c("0" = "Kupffer cells",
               "1" = "Monocytes/cDCs",
               "2" = "Macrophages",
               "3" = "Monocytes/cDCs",
               "4" = "Kupffer cells",
               "5" = "Monocytes/cDCs",
               "6" = "Monocytes/cDCs",
               "7" = "Macrophages",
               "8" = "Monocytes/cDCs",
               "9" = "cDC1",
               "10" = "pDCs",
               "11" = "pDCs",
               "12" = "Dividing cDCs",
               "13" = "Kupffer cells",
               "14" = "Hepatocytes")

all_m_cells$mye_annot = new_m_labs[as.character(all_m_cells$pca25_res.0.6)]

Subset Monocytes

immune_pops = c("2", "7", "1", "6", "8", "3", "5", "4", "13", "0", "9")
all_mon_cells = all_m_cells[,all_m_cells@meta.data$pca25_res.0.6 %in% immune_pops &
                            all_m_cells@meta.data$allcells_clusters %in% c("Macrophages", "cDCs 1",
                                                                           "cDCs 2", "Kupffer cells")]
all_mon_cells = suppressWarnings(SCTransform(all_mon_cells, do.correct.umi = T, verbose = F, 
                                             vars.to.regress=c("unique_name","nCount_RNA"),
                                             variable.features.rv.th = 1, seed.use = 1,
                                             return.only.var.genes = F, 
                                             variable.features.n = NULL))
all_mon_cells = RunPCA(all_mon_cells, verbose = F)
all_mon_cells = RunUMAP(all_mon_cells, dims = 1:25, verbose = F)
DimPlot(all_mon_cells, reduction = "umap", group.by = "Condition")
DimPlot(all_mon_cells, reduction = "umap", group.by = "allcells_clusters")
DimPlot(all_mon_cells, reduction = "umap", group.by = "pca25_res.0.6")
DimPlot(all_mon_cells, reduction = "umap", group.by = "Donor")
DimPlot(all_mon_cells, reduction = "umap", group.by = "Phase")
FeaturePlot(all_mon_cells, reduction = "umap", features = c("MKI67", "ALB", "S100A8", "COLEC11",
                                                            "NKG7", "TRGC1", "PTPRC", "EPCAM",
                                                            "CD14", "TRAC", "CD3E", "CD8A"))

Clustering and markers

immune_pops = c("2", "7", "1", "6", "8", "3", "5", "4", "13", "0", "9")
all_mon_cells = all_m_cells[,all_m_cells@meta.data$pca25_res.0.6 %in% immune_pops &
                            all_m_cells@meta.data$allcells_clusters %in% c("Macrophages", "cDCs 1",
                                                                           "cDCs 2", "Kupffer cells")]
all_mon_cells = suppressWarnings(SCTransform(all_mon_cells, do.correct.umi = T, verbose = F, 
                                             vars.to.regress=c("unique_name","nCount_RNA"),
                                             variable.features.rv.th = 1, seed.use = 1,
                                             return.only.var.genes = F, 
                                             variable.features.n = NULL))

all_mon_cells = RunPCA(all_mon_cells, verbose = F)
all_mon_cells = RunUMAP(all_mon_cells, dims = 1:25, verbose = F)
DimPlot(all_mon_cells, reduction = "umap", group.by = "Condition")

DimPlot(all_mon_cells, reduction = "umap", group.by = "allcells_clusters")

DimPlot(all_mon_cells, reduction = "umap", group.by = "pca25_res.0.6")

DimPlot(all_mon_cells, reduction = "umap", group.by = "Donor")

DimPlot(all_mon_cells, reduction = "umap", group.by = "Phase")

FeaturePlot(all_mon_cells, reduction = "umap", features = c("MKI67", "ALB", "S100A8", "COLEC11",
                                                            "NKG7", "TRGC1", "PTPRC", "EPCAM",
                                                            "CD14", "TRAC", "CD3E", "CD8A"))

Add annotations

DimPlot(all_mon_cells, reduction = "umap", group.by = "pca25_res.0.5", label = T)
DimPlot(all_mon_cells, reduction = "umap", group.by = "allcells_clusters", label = T)

DimPlot(all_mon_cells, reduction = "umap", group.by = "Donor", label = F)

DimPlot(all_mon_cells, reduction = "umap", group.by = "Condition", label = F)

Put annotations on single immune cell object

Make dataframe with new annotations

new_m_labs = c("0" = "Kupffer cells",
               "1" = "cDC2",
               "2" = "Monocytes (IGSF21+ GPR34+)", # similar to those identified here https://www.nature.com/articles/s41586-020-2922-4 
               "3" = "Macrophages (HES4+)", # very similar to 5, HES4 is one of its most unique markers
               "4" = "Kupffer cells (SUCNR1+)", # disproves https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1986575/; these KC are inflammatory/antiviral
               "5" = "Macrophages",
               "6" = "Monocytes (TREM2+ CD9+)", # in ramachandran et al, assoc with fibrotic scars
               ## in their projection, also close to KC
               "7" = "cDC2",
               "8" = "cDC1",
               "9" = "Monocytes (secretory)", # related to 2 and 6
               "10" = "activated DCs") # CD80, CD86, CCR7

all_mon_cells$mono_annot = new_m_labs[as.character(all_mon_cells$pca25_res.0.5)]

Add to the immune Seurat object

all_imm_cells = AddMetaData(all_imm_cells, metadata = newannot_df)
# the original "Dividing cells" will be relabeled "Dividing T/NK cells",
## and the newly annotated "Dividing NK cells" will be renamed to match this
all_imm_cells$immune_annot[is.na(all_imm_cells$immune_annot)] = "Dividing T/NK cells"
all_imm_cells$immune_annot[all_imm_cells$immune_annot=="Dividing NK cells"] = "Dividing T/NK cells"
all_imm_cells$immune_annot[all_imm_cells$immune_annot=="Hepatocytes"] = "Hepatocyte-Monocyte interaction"

DimPlot(all_l_cells, group.by = "lymphoid_annot", reduction = "umap", label = T)+NoLegend()+ggtitle("Lymphoid cells")
DimPlot(all_t_cells, group.by = "t_annot", reduction = "umap", label = T)+NoLegend()+ggtitle("T and NK cells (subset of Lymphoid)")
DimPlot(all_m_cells, group.by = "mye_annot", reduction = "umap", label = T)+NoLegend()+ggtitle("Myeloid cells")
DimPlot(all_mon_cells, group.by = "mono_annot", reduction = "umap", label = T)+NoLegend()+ggtitle("Monocytes cells (subset of Myeloid)")
DimPlot(all_imm_cells, group.by = "allcells_clusters", reduction = "umap", label = T)+NoLegend()+ggtitle("Immune cells (original annotation)")
DimPlot(all_imm_cells, group.by = "immune_annot", reduction = "umap", label = T)+NoLegend()+ggtitle("Immune cells (new detailed annotation)")

Save Seurat objects

all_imm_cells = AddMetaData(all_imm_cells, metadata = newannot_df)
# the original "Dividing cells" will be relabeled "Dividing T/NK cells",
## and the newly annotated "Dividing NK cells" will be renamed to match this
all_imm_cells$immune_annot[is.na(all_imm_cells$immune_annot)] = "Dividing T/NK cells"
all_imm_cells$immune_annot[all_imm_cells$immune_annot=="Dividing NK cells"] = "Dividing T/NK cells"
all_imm_cells$immune_annot[all_imm_cells$immune_annot=="Hepatocytes"] = "Hepatocyte-Monocyte interaction"

DimPlot(all_l_cells, group.by = "lymphoid_annot", reduction = "umap", label = T)+NoLegend()+ggtitle("Lymphoid cells")

DimPlot(all_t_cells, group.by = "t_annot", reduction = "umap", label = T)+NoLegend()+ggtitle("T and NK cells (subset of Lymphoid)")

DimPlot(all_m_cells, group.by = "mye_annot", reduction = "umap", label = T)+NoLegend()+ggtitle("Myeloid cells")

DimPlot(all_mon_cells, group.by = "mono_annot", reduction = "umap", label = T)+NoLegend()+ggtitle("Monocytes cells (subset of Myeloid)")

DimPlot(all_imm_cells, group.by = "allcells_clusters", reduction = "umap", label = T)+NoLegend()+ggtitle("Immune cells (original annotation)")

DimPlot(all_imm_cells, group.by = "immune_annot", reduction = "umap", label = T)+NoLegend()+ggtitle("Immune cells (new detailed annotation)")